Aulas de Introdução a Ciências da Computação (2020-02)¶
Aula 01¶
Análise: Base de dados Iris¶
plot(iris[,1:4], col = iris[,5], pch=20, cex=1.3)
Correlações de variáveis¶
names(iris)
cor(iris[,1], iris[,2])
cor(iris[,1], iris[,3])
Produto escalar entre matrizes:
iris_1_centralizada = iris[,1]-mean(iris[,1])
iris_2_centralizada = iris[,2]-mean(iris[,2])
iris_1_centralizada%*%iris_2_centralizada
x = sin(2*pi*seq(0,9,len=1000))
y = x + rnorm(mean=0, sd=.1, n=1000)
plot(x, cex=1.2, pch=20)
points(y, col=2, cex=.6, pch=20)
Similaridade média entre pontos
((x-mean(x))%*%(y-mean(y)) / (sd(x)*sd(y))) / 1000
cor(x,y)
plot(iris[,3:4], col=iris[,5], cex=1.6, pch=20)
Aula 02¶
Análise de PCA¶
install.packages("splus2R")
require(splus2R)
data = as.data.frame(rmvnorm(100, rho=.9) * 5+70)
colnames(data) = c("Statistics", "Physics")
plot(data, pch=20, cex=1.5)
scaled.cov = cov(scaled)
print(scaled.cov)
eigenspace = eigen(scaled.cov)
print(eigenspace)
rownames(eigenspace$vectors) = c("Statistics", "Physics")
colnames(eigenspace$vectors) = c("PC1", "PC2")
print(eigenspace)
scaled = apply(data, 2, function(x) {x - mean(x)})
plot(scaled, pch=20, cex=1.5)
points(scaled[21], col=2, pch=20, cex=1.5)
x2 = apply(data, 2, function(col) {(col - mean(col)) / sd(col)})
points(x2, col=2, pch=20, cex=1.5)
abline(v=0,h=0)
pc1.slope = eigenspace$vectors[1,1]/eigenspace$vectors[2,1]
pc2.slope = eigenspace$vectors[1,2]/eigenspace$vectors[2,2]
abline(0, pc1.slope, col='green')
abline(0, pc2.slope, col='blue')
sum(eigenspace$values)
sum(var(scaled[,1]), var(scaled[,2]))
eigenspace$values/sum(eigenspace$values)
Score plot¶
scores = scaled %*% eigenspace$vectors
plot(scores)
abline(0,0, col='red')
abline(0,90, col='blue')
par(mfrow = c(2,1))
plot(scaled, pch=20, cex=1.5)
abline(0, pc1.slope, col='red')
abline(0, pc2.slope, col='blue')
plot(scores, pch=20, cex=1.5)
abline(0,0, col='red')
abline(0,90, col='blue')
Biplot¶
Com essa forma de gráfico, podemos entender melhor a relação existente entre as variáveis originais
plot(scores, xlim=range(scores), ylim=range(scores), col='orange')
abline(0,0, col='red')
abline(0,90, col='blue')
sd = sqrt(eigenspace$values)
factor = 1
arrows(0,0, eigenspace$vectors[,1] * sd[1]/factor, eigenspace$vectors[,2]*sd[2]/factor, length=.1, col=1)
text(eigenspace$vectors[,1]*sd[1]/factor*1.2,
eigenspace$vectors[,2]*sd[2]/factor*1.2, c("Statistics", "Physics"), cex=.8, col=2)
Como o ângulo entre os vetores é agudo, existe forte correlação entre as variáveis, e pontos entre esses vetores apresentam um equilíbrio entre cada variável.
De outra forma, pontos acima da linha vermelha, cujo ângulo com o vetor "Statistics" é menor do que com o vetor "Physics", são pontos que apresentam maiores notas em estatística.
Análise: Base de dados Wine¶
wine_data = read.table("Estudos/WineData/wine.data", sep=',')
wine_data[0:3,]
Y = wine_data[,1]
X = wine_data[,2:14]
plot(X[,5:8], col=Y, pch=20, cex=2)
Variação de escala entre os parâmetros observados (colunas)
apply(X, 2, function(col) {range(col)})
Aplicando correção de centralização e reescala
scaled = apply(X, 2, function(col) {(col - mean(col)) / sd(col)})
apply(scaled, 2, function(col) {range(col)})
scaled[2:5,]
scaled.cov = cov(scaled)
eigenspace = eigen(scaled.cov)
eigenspace
Onde cada coluna representa uma componente principal, e cada linha, o vetor resultante de cada atributo
eigenspace$values / sum(eigenspace$values)
plot(cumsum(eigenspace$values / sum(eigenspace$values)), pch=20, cex=1.5)
abline(v=8, h=.9)
print(as.numeric(scaled[1,]))
Exemplo: aplicando o produto escalar sobre os autovetores e a primeira coluna do conjunto de dados escalados
print((eigenspace$vectors)%*%(as.numeric(scaled[1,])))
Aplicando o produto escalar sobre todo o conjunto de dados
scores = scaled %*% eigenspace$vectors
scores = as.data.frame(scaled %*% eigenspace$vectors)
colnames(scores) = paste("PC", 1:13, sep='')
print(colnames(scores))
plot(scores[,1:4], col=Y, pch=20, cex=1.5)
Análise KNN: K-Nearest Neighbours¶
Considerando a média do atribudo dos $k$ vizinhos mais próximos
knn <- function(query, k, X, Y) {
#Distancia euclidiana que a consulta (query) tem em relação aos exemplos (linhas) da matriz X
E = apply(X, 1, function(row) { sqrt(sum((row-query))^2) })
#Posição relativa dos k vizinhos mais próximos
ids = sort.list(E, dec=F)[1:k]
# Classes únicas
classes = unique(Y)
i = 1
count = rep(0, length(classes))
for (class in classes) {
count[i] = sum(class == Y[ids])
i = i + 1
}
ret = list()
ret$classes = classes
ret$count = count
ret$max.prob.class = classes[which.max(count)]
return (ret)
}
wine_data = read.table("Estudos/WineData/wine.data", sep=',')
Y = wine_data[,1]
X = wine_data[,2:14]
res = knn(query=X[1,], k=5, X=X, Y=Y)
print(res)
res = knn(query=X[1,], k=50, X=X, Y=Y)
print(res)
Conjunto de treinamento
ids = sample(1:nrow(X), size=ceiling(0.5*nrow(X)))
X.train = X[ids,]
Y.train = Y[ids]
print(X.train[,3])
print("Classes:")
print(Y.train[0:3])
Conjunto de teste
X.test = X[-ids,]
Y.test = Y[-ids]
# X[-ids] <-- entende que -ids representa o número das colunas (e pega todas as linhas)
# X[-ids,] <-- entende que -ids representa o número das linhas (e pega todas as colunas)
accuracy = 0
for (i in 1:nrow(X.test)) {
res = knn(query=X.test[i,], k=3, X=X.train, Y=Y.train)
res_class = res$max.prob.class
y_class = Y.test[i]
if (res_class == y_class) {
accuracy = accuracy + 1
}
}
accuracy = accuracy / nrow(X.test)
cat("Accuracy = ", accuracy, "\n")
R = NULL
for (k in 1:10) {
accuracy = 0
for (i in 1:nrow(X.test)) {
res = knn(query=X.test[i,], k=k, X=X.train, Y=Y.train)
res_class = res$max.prob.class
y_class = Y.test[i]
if (res_class == y_class) {
accuracy = accuracy + 1
}
}
accuracy = accuracy / nrow(X.test)
R = rbind(R, c(k, accuracy))
}
plot(R, xlab='k values', ylab='Accuracy', pch=20, cex=1.5)
lines(R, pch=20, cex=1.5)
Aula 03¶
Aprendizado Baseado em Instâncias¶
Combinando os resultados anteriores em uma única função
test.wine <- function(PCA = FALSE, perc = 0.9) {
# PCA == FALSE --> Sem transformação por auto-valores/auto-vetores
# PCA == TRUE --> Com transformação por PCA em um subconjunto de PCs
wine_data = read.table("Estudos/WineData/wine.data", sep=',')
Y = wine_data[,1]
X = wine_data[,2:14]
if (PCA) {
# Transformação espacial
# Centralização e reescala
scaled = apply(X, 2, function(col) {(col - mean(col)) / sd(col)})
# Avaliação das dependências lineares entre os atributos
C = cov(scaled)
E = eigen(C)
# Plot da soma acumulada de autovalores
plot(cumsum(E$values/sum(E$values)), main="Importância agregada por PC", xlab='# PCs',
ylab='Importância em %',, pch=20, cex=1.5)
# Returnar o número de autovalores tais que a soma acumulada seja maior do que um dado percentual
numero.PCs = which(cumsum(E$values/sum(E$values)) >= perc)[1]
print(numero.PCs)
# Realizar transformação do espaço (projeção dos dados sobre os autovetores)
Transform = scaled %*% E$vectors
X = Transform[,1:numero.PCs]
}
R = matrix(0, nrow=10, ncol=5)
for (j in 1:5){
#Hold-out
# - uma parcela dos dados é escolhida aleatoriamente (uniforme) para
# para compor o conjunto de treinamento e o restante disjunto dos exemplos irá
# compor o conjunto de treinamento
ids = sample(1:nrow(X), size=ceiling(0.75*nrow(X)))
X.train = X[ids,]
Y.train = Y[ids]
X.test = X[-ids,]
Y.test = Y[-ids]
for (k in 1:10) {
accuracy = 0
for (i in 1:nrow(X.test)) {
res = knn(query=X.test[i,], k=k, X=X.train, Y=Y.train)
res_class = res$max.prob.class
y_class = Y.test[i]
if (res_class == y_class) {
accuracy = accuracy + 1
}
}
accuracy = accuracy / nrow(X.test)
R[k,j] = accuracy
}
}
return (R)
}
com.PCA = test.wine(PCA=TRUE, perc=.95)
Analisando a média dos valores de acurácia obtidos com Hold-out
apply(com.PCA, 1, function(row) {mean(row)})
Análise DWNN: Distance-Weighted Nearest Neighbors¶
KNN -> Sistema de contagem (classes discretas) para problemas de categorização
DWNN -> Regressão do espaço (classes contínuas) para problemas de regressão
No caso discreto, tomamos o argumento máximo (análogo a uma heaviside)
$$ \hat{f}(X_q) = argmax_{v \in V} \sum_{i=1}^{k}w_i \delta(v, f(x_i) $$$$ w_i \equiv \frac{1}{d(x_q, x_i)^2} $$com $d(a,b_i)$ distância euclidiana entre o ponto de consulta $a$ e os pontos vizinhos $b_i$
No caso contínuo
$$ \hat{f}(x_q) \leftarrow \frac{\sum_{i=1}^{k} w_i f(x_i)}{\sum_{i=1}^k w_i} $$$$ w_i \equiv \frac{1}{d(x_q, x_i)^2} $$A classe do ponto de consulta é dada pela média ponderada das saídas produzidas pelos vizinhos mais próximos, mapeados por gaussianas centradas no ponto de consulta $x_q$
Neste caso, consideramos um \textit{kernel} radiam tal que os pesos $w_i$ são dados por
$$ w_i = e^{-\frac{E(x_q, x_i)^2}{2 \sigma ^2}} $$ps = c(1,10)
sigma = 1
exp(-ps^2 / (2*sigma))
Como o centro da gaussiana (média) está em zero, a distância até o primeiro ponto de consulta (=1) apresenta importância de 0.6, entquanto um ponto com distância 10, é zero (erro numérico). Caso sigma aumente, a variância da gaussiana no ponto de consulta aumenta e as distâncias são consideradas como mais importantes (a abertura da gaussiana é maior)
sigma = 10
exp(-ps^2 / (2*sigma))
Considerando os pesos para os pontos 1 e 10 como: w(1) = 0.5, w(10) = 127,87
sigma = 1
w = exp(-ps^2 / (2*sigma))
y = c(0.5, 127.87)
w*y
Temos a importância ponderada por cada ponto de acordo com sua distância do ponto de consulta no centro da gaussiana
sum(w*y) / sum(w)
E a soma retorna uma regressão numérica, baseada no peso de cada pontoe em sua distância relativa
dwnn <- function(query, X, Y, sigma) {
E = apply(X, 1, function(row) { sqrt(sum((row - query)^2)) })
# Construção de gaussiana centrada na query
weight = exp(-E^2 / (2*sigma^2))
# Multiplicando ponto a ponto os pesos pela saída de cada elemento
return (weight %*% Y / sum(weight))
}
test.linear <- function(sigma) {
pdf("dwnn_linear_eg.pdf", onefile = FALSE)
dataset = cbind(-5:5, -5:5)
print(dataset)
plot(dataset, pch=20, cex=2, cex.lab=2, cex.axis=2)
X = matrix(dataset[,1], ncol=1)
Y = matrix(dataset[,2], ncol=1)
R = NULL
for (x in seq(-5,5,len=100)) {
y = dwnn(query=x, X=X, Y=Y, sigma=sigma)
R = rbind(R, c(x,y))
}
points(R, col=2, pch=20)
dev.off()
}
Caso $\sigma$ -> 0, apenas pontos pertencentes ao conjunto de dados pode ser analisado, já que o raio da gaussiana é pequeno o suficiente para que não exista vizinhos próximos
Caso $\sigma$ -> $\infty$, a resposta da regressão retorna a média de todos os pontos, já que a gaussiana ocupa todo o espaço e todos os pontos são considerados vizinhos
test.linear(0.002)
test.linear(0.2)
test.linear(0.5)
test.linear(10)
Quanto maior a abertura da gaussiana centrada em zero, todos os pontos do espaço serão considerados com o mesmo peso, mesmo que estejam a -5 e +5 (extremos) de distância do ponto de consulta
test.linear2 <- function(sigma) {
dataset = cbind(-5:5, -5:5)
print(dataset)
plot(dataset, pch=20, cex=2, cex.lab=2, cex.axis=2, xlim=c(-7,7))
X = matrix(dataset[,1], ncol=1)
Y = matrix(dataset[,2], ncol=1)
R = NULL
for (x in seq(-7, 7,len=100)) {
y = dwnn(query=x, X=X, Y=Y, sigma=sigma)
R = rbind(R, c(x,y))
}
points(R, col=2, pch=20)
}
test.linear2(0.2)
Nesse caso, a distribuição de probabilidades conjunta P(x,y) no treinamento (pontos pretos) é distinta da distribuição do teste (pontos vermelhos) passando dos extremos
Regressão para uma série temporal¶
test.sin <- function(sigma) {
# Como serie temporal, é uma variável unidimensional
time = seq(0,9, length=1000)
series = sin(sin(2*pi*time))
par(mfrow=c(2,1))
plot(series, xlim=c(0,1200), pch=20, cex=1, cex.lab=2)
# Temos que adicionar um espaço de entradas X e um espaço de saídas Y
# Dado series[1], esperamos ter o resultado series[2]
X = matrix(series[1:999], ncol=1)
Y = as.numeric(series[2:1000])
# Conjunto de treinamento
ids = 1:800
X.train = matrix(X[ids,], ncol=1)
Y.train = Y[ids]
# Conjunto de teste
X.test = matrix(X[-ids,], ncol=1)
Y.test = Y[-ids]
# Teste
R = NULL
for (i in 1:nrow(X.test)) {
x = X.test[i,]
y = dwnn(query=x, X=X.train, Y=Y.train, sigma=sigma)
R = rbind(R, c(x,y))
}
plot(series[-ids], pch=20, cex=1, cex.lab=2)
points(R[,2], pch=20, cex=.6, col=2)
}
test.sin(2)
test.sin(.6)
test.sin2 <- function(sigma, lag.max=200) {
# Como serie temporal, é uma variável unidimensional
time = seq(0,9, length=1000)
series = sin(sin(2*pi*time))
plot(series, xlim=c(0,1200), pch=20, cex=1, cex.lab=2)
# Temos que adicionar um espaço de entradas X e um espaço de saídas Y
# Dado series[1], esperamos ter o resultado series[2]
X = matrix(series[1:999], ncol=1)
Y = as.numeric(series[2:1000])
# Conjunto de treinamento
X.train = X
Y.train = Y
# Predição recorrente ---> Redes neurais profundas
# - Uma saída vira entrada para prever umapróxima saída
R = matrix(NA, nrow=nrow(X.train), ncol=2)
x = series[1000]
for (i in 1:lag.max) {
y = dwnn(query=x, X=X.train, Y=Y.train, sigma=sigma)
R = rbind(R, c(x,y))
x = y
}
points(R[,2], pch=20, cex=.6, col=2)
}
test.sin2(200)
Caso a abertura seja muito grande, como anteriormente, os dados preditos tendem à média de todos os pontos
test.sin2(.0001)
Considerando redução da abertura da gaussiana, a predição não tende à média.
A distribuição de probabilidades conjunta P(X,Y) para a função test.linear não para a test.sin2
Aula 04¶
Retomada dos algoritmos KNN e DWNN - Paradigmas do IBL (Instance-Based Learning)
knn <- function(query, k, X, Y) {
E = apply(X, 1, function(row) { sqrt(sum((row - query)^2)) })
ids = sort.list(E, dec=F)[1:k]
classes = unique(Y)
count = rep(0, length(classes))
i = 1
for (class in classes) {
count[i] = sum(class == Y[ids])
i = i + 1
}
ret = list()
ret$classes = classes
ret$count = count
ret$max.prob.class = classes[which.max(count)]
return (ret)
}
dwnn <- function(query, X, Y, sigma) {
E = apply(X, 1, function(row) { sqrt(sum((row - query)^2)) })
weight = exp(-E^2 / (2*sigma^2))
return (weight %*% Y / sum(weight))
}
id = c(sample(1:50, size=25), sample(51:100, size=25), sample(101:150, size=25))
print(iris[125,])
query = iris[125,1:4]
knn(query=query, k=15, X=iris[id,1:4], Y=iris[id,5])
Exemplo de uso: DWNN Sunspots
series = read.table("Estudos/SunspotsTimeSeries/monthly-sunspots.csv", sep=',')
series[0:5,]
install.packages("tseriesChaos")
require(tseriesChaos)
dataset = tseriesChaos::embedd(series, m=3, d=5)
plot(c(1:length(series[,1])), series[,2])
Esse pacote permite a criação de um dataset que separa os dados em três colunas (d=3) relativas a um delay temporal de 5 unidades (m=5). Assim, podemos usar o dataset e a função dwnn para prever dados ao longo do tempo
dataset[1:10,]
dataset[2001,1:2]
dwnn(query=dataset[2001,1:2], X=dataset[1:2000, 1:2], Y=dataset[1:2000, 3], sigma=1)
dwnn(query=dataset[2001,1:2], X=dataset[1:2000, 1:2], Y=dataset[1:2000, 3], sigma=5)
dwnn(query=dataset[2001,1:2], X=dataset[1:2000, 1:2], Y=dataset[1:2000, 3], sigma=50)
Análise com KNN¶
Dividir o conjunto em treinamento e teste:
- Hold-out: separar o conjunto original em 2 subconjunto disjuntos: treinamento e teste (+grosseiro)
- Leave-one-out: separar um único elemento do original como teste. O resto é treinamento (+preciso)
- K-fold Cross Validation: dividir o conjunto original em $k$ subconjuntos. Treinamento com $k-1$. Teste com $k-ésimo$ (balanceado)
Visualizaros resultados de classificação
Aula 05¶
dwnn <- function(query, X, Y, sigma) {
E = apply(X, 1, function(row) { sqrt(sum((row - query)^2)) })
# Construção de gaussiana centrada na query
weight = exp(-E^2 / (2*sigma^2))
# Multiplicando ponto a ponto os pesos pela saída de cada elemento
return (weight %*% Y / sum(weight))
}
Mesmo que a função não seja definida para valores intermediários entre os pontos inteiros de [0,10], com o DWNN podemos tomar esses valores como pontos de consult, tanto para teste quanto para treinamento
test.linear.tempo <- function(sigma=0.5) {
# Tempo Valor da observação
dataset = cbind(seq(0,10,length=50), seq(20,25, length=50))
plot(dataset, cex=1.6, pch=20, cex.axis=2, xlab='Tempo', ylab='Observável', xlim=c(-5,15))
X = matrix(dataset[,1], ncol=1)
Y = dataset[,2]
N = NULL
for (tempo in seq(-5, 20, length=150)) {
valor = dwnn(query=tempo, X=X, Y=Y, sigma=sigma)
N = rbind(N, c(tempo, valor))
}
points(N, col=2, pch=20, cex=.8)
}
test.linear.tempo()
Nesse sentido, o espaço de valores possíveis é considerado como apenas aqueles pertencentes ao conjunto de pontos existentes na minha amostra (nesse caso, de 0 a 10). Logo a extrapolação de valores além desse limite gera predição igual ao valor extremo do intervalo
Sendo $P(X, Y)$ probabilidade conjunta entre as variáveis $X$ e $Y$
$P(X,Y)$ -> contém a probabilidade associada a cada evento discreto. Então qualquer ponto que vá além do espaço daamostra não consegue ser previsto
test.seno.tempo <- function(sigma=0.5) {
tempo = seq(0, 3, length=250)
dataset = cbind(tempo, sin(2*pi*tempo))
plot(dataset, cex=1.6, pch=20, cex.axis=2, xlab='Tempo', ylab='Observável', xlim=c(0,4))
X = matrix(dataset[,1], ncol=1)
Y = dataset[,2]
N = NULL
for (tempo in seq(0, 4, length=350)) {
valor = dwnn(query=tempo, X=X, Y=Y, sigma=sigma)
N = rbind(N, c(tempo, valor))
}
points(N, col=2, pch=20, cex=.8)
}
test.seno.tempo()
test.seno.tempo(0.1)
Reorganizando o conjunto de dados em uma outra forma de espaço
Takens embedding theorem¶
Há um espaço que permite mapear uma observação x(t) em um determinado tempo a demais observações x(t-d), x(t-2d), ... supondo que o sistema seja determinístico
tempo = seq(0, 3, length=250)
valores = sin(2*pi*tempo)
train = valores[1:125]
test = valores[126:250]
espaco.train = cbind(train[1:124], train[2:125])
espaco.test = cbind(test[1:124], test[2:125])
espaco.train[1:5,]
plot(espaco.train, xlab='x(t)', ylab='x(t+1)', cex=1.4, pch=20)
points(espaco.test, c=2 , cex=.4, pch=20)
Agora a distribuição de probabilidade conjunta $P(X,Y)$, no espaço reconstruído considera as relações que ocorrem em instantes de tempos distintos
Caso seja determinístico, os dados para prever o futuro são condizentes com os dados do passado, e dado um valor de teste x(t), a predição para x(t+1) é coerente
require(tseriesChaos)
test.seno.espaco_fase <- function(sigma=0.5) {
tempo = seq(0, 3, length=250)
calores = sin(2*pi*tempo)
dataset = tseriesChaos::embedd(valores, m=2, d=1)
X = matrix(dataset[,1], ncol=1)
Y = dataset[,2]
x.last = as.numeric(Y[length(Y)])
N = NULL
for (i in 1:50) {
x.last.plus.1 = dwnn(query=x.last, X=X, Y=Y, sigma=sigma)
N = rbind(N, c(x.last, x.last.plus.1))
x.last = x.last.plus.1
}
plot(dataset, cex=1.4, pch=20, cex.lab=2, cex.axis=1.5, xlab='x(t)', ylab='x(t+1)')
points(N, col=2, pch=20, cex=.8)
}
test.seno.espaco_fase()
test.seno.espaco_fase(.01)
test.seno.espaco_fase(.001)
Série temporal: sequência de valores obtidos ao longo do tempo (uniformemente distribuídos ou não)
- Puramente estocástica (sem dependência temporal)
- Puramente determinística (todo valor pode ser representado por uma combinação de valores históricos)
- Misturas
Analisando a correlação da função sobre si mesma
series = sin(2*pi*seq(0,9,len=1000))
series[1:1000] %*% series[1:1000]
series[1:999] %*% series[2:1000]
series[1:998] %*% series[3:1000]
acf(series, lag.max=100)
Para cada valor ao longo do eixo "Lag", "ACF" (Auto-Correlation Function) representa a correlação da função sobre si mesma
Para valores fora do intervalo de confiança, existe dependência temporal significativa
Caso haja dependência, oTeorema de Takens pode ser aplicado para modelar o relação temporal
Utilizar um algoritmo de regressão (DWNN, e.g.)
Caso resultados ruins, decompor em componentes estocástica e determinística
Artigo recomendado: https://journals.aps.org/pra/abstract/10.1103/PhysRevA.33.1134
mutual(series, lag.max = 50)
Os mínimos da informação mútua são importantes candidatos a o número de lag temporal que otimiza a predição via Takens embedding
Artigo recomendado: https://journals.aps.org/pra/abstract/10.1103/PhysRevA.45.3403
plot(false.nearest(series, d=5, m=10, t=10))
plot(false.nearest(series, d=16, m=10, t=10))
test.seno_com_ruido.espaco_fase <- function(sigma=0.5) {
tempo = seq(0, 5, length=1000)
valores = sin(2*pi*tempo) + rnorm(mean=0, sd=0.15, n=1000)
par(mfrow=c(2,1))
plot(valores)
acf(valores)
}
test.seno_com_ruido.espaco_fase()
Pelo resultado do ACF: Tem dependência temporal expressiva
valores[0:10]
z = mutual(valores[0:10])
print(z)
print(diff(z))
false.nearest(valores, m=10, d=5, t=10)
print(false.nearest(valores, m=10, d=5, t=10)[1,])
test.seno_com_ruido.espaco_fase <- function(sigma=0.5, sd=0.1) {
tempo = seq(0, 5, length=1000)
valores = sin(2*pi*tempo) + rnorm(mean=0, sd=sd, n=1000)
# Tentativa de modelagem determinística
ami = tseriesChaos::mutual(valores, lag.max=100)
# ami <--- 1 0.8 0.7 0.6 0.7
v = diff(ami)
# v <----- -0.2 -0.1 -0.1 +0.1
d = as.numeric(which(v > 0)[1]) # Primeiro elemento cujo valor ami é positivo
fnn = tseriesChaos::false.nearest(valores, m=10, d=d, t=10)
print(fnn)
m = as.numeric(which.min(fnn[1,])[1])
cat("Embedding dimension = ", m, "\n")
cat("Time lag = ", d, "\n")
dataset = tseriesChaos::embedd(valores, m=m, d=d)
labelId = ncol(dataset)
print(labelId)
X = matrix(dataset[,1:(labelId-1)], ncol=labelId-1)
Y = dataset[,labelId]
x.last = as.numeric(Y[length(Y)])
N = NULL
buffer = dataset
for (new.row in (nrow(buffer)+1):(nrow(buffer)+300)){
x = as.numeric(buffer[new.row - d, 2:ncol(buffer)])
y = dwnn(query=x, X=X, Y=Y, sigma=sigma)
buffer = rbind(buffer, c(x,y))
}
plot(Y, t='l', cex=1.5, pch=20, xlab='Tempo', ylab='Valores', xlim=c(1,nrow(buffer)))
points(buffer[,labelId], col=2, pch=20, cex=.5)
}
test.seno_com_ruido.espaco_fase(0.2)
test.seno_com_ruido.espaco_fase(0.2, 0.2)
test.seno_com_ruido.espaco_fase(0.08, 0.2)
Recomendação de artigo: https://www.sciencedirect.com/science/article/abs/pii/S0165168415002297?via%3Dihub
Aula 06¶
dwnn <- function(query, X, Y, sigma) {
E = apply(X, 1, function(row) { sqrt(sum((row - query)^2)) })
# Construção de gaussiana centrada na query
weight = exp(-E^2 / (2*sigma^2))
# Multiplicando ponto a ponto os pesos pela saída de cada elemento
return (weight %*% Y / sum(weight))
}
test.seno_com_muito_ruido.espaco_fase <- function(sigma=0.5, sd=0.5) {
tempo = seq(0,5, length=1000)
valores = sin(2*pi*tempo) + rnorm(mean=0, sd=sd, n=1000)
par(mfrow=c(3,1))
ami = tseriesChaos::mutual(valores, lag.max=100)
v = diff(ami)
d = as.numeric(which(v > 0)[1])
cat("Time lag = ", d, "\n")
# Estimando Takens embedding dimension
fnn = tseriesChaos::false.nearest(valores, m=10, d=d, t=10)
plot(fnn)
print(fnn)
m = as.numeric(which.min(fnn[1,]))
# Aplicação do Takens embedding theorem
dataset = tseriesChaos::embedd(valores, m=m, d=d)
labelId = ncol(dataset)
X = matrix(dataset[,1:(labelId-1)], ncol=labelId-1)
Y = dataset[,labelId]
buffer = dataset
for (new.row in (nrow(buffer) + 1) :(nrow(buffer) + 300)) {
x = as.numeric(buffer[new.row - d, 2:ncol(buffer)])
y = dwnn(query=x, X=X, Y=Y, sigma=sigma)
buffer = rbind(buffer, c(x,y))
}
plot(Y, t='l', cex.lab=2, cex=2, cex.axis=2,
xlab="Tempo", ylab='Valores', xlim=c(1, nrow(buffer)))
points(buffer[, labelId], col=2, pch=20, cex=.2)
}
test.seno_com_muito_ruido.espaco_fase(sigma=0.5, sd=0.05)
O resultado precisa ter no mínimo percentual abaixo de 20% (como mencionado no artigo que deu origem ao método de False Nearest Neighbors)
test.seno_com_muito_ruido.espaco_fase(sigma=0.0005, sd=0.5)
Neste cenário com muito ruído, existe um percentual alto suficiente de falsos vizinhos (acima de 70%) que o modelo nem ao menos consegue prever e não encontra vizinho algum. Neste sentido, a solução que temos é através da decomposição da senóide em componentes determinístico e estocástico
Consideramos a adição entre componente determinístico e estocástico
Estocástico: $e(t) = N(\mu = 0, \sigma = 0.5)$
Determinístico: $d(t) = sin(2\pi t)$
Componente final: $x(t) = d(t) + e(t)$
install.packages("EMD")
require(EMD)
EMD: Empirical Mode Decomposition -> Criação de um envelope tomando os máximos e mínimos de uma dada série temporal. Em seguida, subtrai a função superior e inferior e toma, assim, uma função média para os dados. Essa função média pode ser então comparada com o valor real do dado naquele ponto
test.seno_com_muito_ruido.decomposicao.espaco_fase <- function(sigma=0.5, sd=0.5) {
tempo = seq(0,5, length=1000)
valores = sin(2*pi*tempo) + rnorm(mean=0, sd=sd, n=1000)
# Componente estocástico
# Componente determinístico
return (valores)
}
x = test.seno_com_muito_ruido.decomposicao.espaco_fase()
plot(x, pch=20, cex=.6)
media = 0
m_t = c()
# Usando média móvel ponderada por função exponencial
for (i in 1:length(x)) { media = media * 0.9 + x[i]*0.1; m_t=c(m_t, media); }
lines(m_t, col=2)
new_x = x - m_t
lines(new_x, col=3)
A nova função (subtração do dado real com a média móvel ponderada), em verde, permite analisar:
Zero-crossing: Se o número de vezes que a função atravessou a linha de zero é pequeno, significa que não existe uma tendência de senóide
Isso permite decompora função original em senóides
decomposicao = emd(x)
decomposicao$nimf # Number Intrisic Mode Function (Número de senóides que é encontrada na série original)
par(mfrow=c(3,3))
plot(x, cex.axis=1, t='l')
plot(decomposicao$imf[,1], cex.axis=1, t='l') # Primeira sinusoidal encontrada
plot(decomposicao$imf[,2], cex.axis=1, t='l') # Segunda sinusoidal encontrada
plot(decomposicao$imf[,3], cex.axis=1, t='l')
plot(decomposicao$imf[,4], cex.axis=1, t='l')
plot(decomposicao$imf[,5], cex.axis=1, t='l')
plot(decomposicao$imf[,6], cex.axis=1, t='l')
plot(decomposicao$imf[,7], cex.axis=1, t='l')
plot(decomposicao$residue, cex.axis=1, t='l')
O resíduo é o conjunto decomposto que não conseguiu completar um período de senóide
par(mfrow=c(2,1)); plot(decomposicao$imf[,1], cex.axis=1, pch=20)
plot(fft(decomposicao$imf[,1]), asp=1, cex=.8, pch=20) # Transformada de Fourier para a primeira decomposição
par(mfrow=c(3,3))
plot(fft(decomposicao$imf[,1]), asp=1, cex=.8, pch=20) # Transformada de Fourier para a primeira decomposição
plot(fft(decomposicao$imf[,2]), asp=1, cex=.8, pch=20) # Transformada de Fourier para a segunda decomposição
plot(fft(decomposicao$imf[,3]), asp=1, cex=.8, pch=20)
plot(fft(decomposicao$imf[,4]), asp=1, cex=.8, pch=20)
plot(fft(decomposicao$imf[,5]), asp=1, cex=.8, pch=20)
plot(fft(decomposicao$imf[,6]), asp=1, cex=.8, pch=20)
plot(fft(decomposicao$imf[,7]), asp=1, cex=.8, pch=20)
plot(fft(decomposicao$residue), asp=1, cex=.8, pch=20)
Quando tivermos vetores da decomposição de Fourier que apresentem norma muito grande (Caso da sexta componente EMD), é um bom sinal de que começamos a decomposição determinística.
Além disso: Congruência de Fase: Quando um vetor apresenta vizinhos próximos a ele que apresentam período ao redor do plano complexo onde apenas o angulo desse vetor altera de forma mais suave.
Sendo coeff = atan(Im(coeff) / Re(coeff)
o ângulo entre as componentes Real e Imaginária do plano complexo, podemos:
coeff = fft(decomposicao$imf[,1])
par(mfrow=c(3,1))
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20) # Mapa dos ângulos da primeira decomposição em senóide
coeff = fft(decomposicao$imf[,2])
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20)
coeff = fft(decomposicao$imf[,3])
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20)
Começa a aparecer uma congruência de fase ao centro. São angulos existentes entre os vetores do espaço de Fourier que começama transicionar de maneira mais suave
coeff = fft(decomposicao$imf[,1])
par(mfrow=c(3,3))
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20) # Mapa dos ângulos da primeira decomposição em senóide
coeff = fft(decomposicao$imf[,2])
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20)
coeff = fft(decomposicao$imf[,3])
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20)
coeff = fft(decomposicao$imf[,4])
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20)
coeff = fft(decomposicao$imf[,5])
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20)
coeff = fft(decomposicao$imf[,6])
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20)
coeff = fft(decomposicao$imf[,7])
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20)
coeff = fft(decomposicao$residue)
plot(atan(Im(coeff) / Re(coeff)), cex=.8, pch=20)
Uma das conclusões: A primeira decomposição apresenta baixa congruência de fase, que vai aumentando conforme as decomposições avançam
estocastico = rowSums(decomposicao$imf[, 1:4])
deterministico = rowSums(decomposicao$imf[, 5:7]) + decomposicao$residue
par(mfrow=c(3,1))
plot(x, pch=20, cex=.8)
plot(estocastico, pch=20, cex=.8)
plot(deterministico, pch=20, cex=.8)
Com essa decomposição, podemos modelar de forma mais assertiva o a série temporal dada. A componente estocástica pode ser modelada por ferramental estocástico e determinístico, com ferramental de processos dinâmicos
t.test(estocastico, rnorm(mean=0, sd=0.5, n=1000)) # Usando de teste de hipótese para
# comparar a nossa função original
O p-value pode ser também compreendido como a área de interseção entre dois histogramas representados pelas funções que queremos comparar pelo test t
par(mfrow=c(2,1))
hist(estocastico, breaks=50)
hist(rnorm(mean=0, sd=0.5, n=1000), breaks=50)
Assim, um valor de p-value = 0.758 representa que 75.8% das áreas dos histogramas se interceptam
Agora precisamos comparar uma congruência de fase com as anteriores e saber o quanto de ganho tivemos. Podemos fazer isso através da Mutual Information (Diferente do que já vimos, que se refere a AUTO Mutual Information)
install.packages('c3net')
require(c3net)
makemim(rbind(rnorm(mean=0, sd=1, n=1000), rnorm(mean=0, sd=1, n=1000)))
makemim(rbind(sin(2*pi*seq(0,9,len=1000)), sin(2*pi*seq(0,9,len=1000))))
A informação mútua entre duas senóides é alta.
A mutual information considera a entropia de duas variáveis aleatórias, sendo assim a entropia conjunta de $X$ e $Y$, mostrando qual a dependência (não necessariamente linear, que é o caso auto-mutual information) entre duas variáveis
test.seno_com_muito_ruido.decomposicao.espaco_fase <- function(sigma=0.5, sd=0.5, mi.epsilon=0.5) {
tempo = seq(0,5, length=1000)
valores = sin(2*pi*tempo) + rnorm(mean=0, sd=sd, n=1000)
# Componente estocástico
require(EMD)
decomposicao = emd(valores)
phases.vector = list()
for (i in 1:decomposicao$nimf) {
coeff = fft(decomposicao$imf[,i])
phases = atan(Im(coeff) / Re(coeff))
phases.vector[[i]] = phases
}
coeff = fft(decomposicao$residue)
phases = atan(Im(coeff) / Re(coeff))
phases.vector[[i+1]] = phases
require(c3net)
mutual.information = rep(0, length(phases.vector)-1)
for (i in 1:(length(phases.vector) - 1)) {
mutual.information[i] = makemim(rbind(phases.vector[[i]], phases.vector[[i+1]]))[1,2]
}
plot(mutual.information, pch=20, cex=.8)
det.start = which(mutual.information > mi.epsilon)[1]
stochastic = rowSums(decomposicao$imf[,1:(det.start-1)])
# Componente determinístico
deterministic = rowSums(decomposicao$imf[,det.start:decomposicao$nimf]) + decomposicao$residue
par(mfrow=c(3,1))
plot(valores, pch=20, cex=.8, cex.axis=2, main='Serie temporal original')
plot(stochastic, pch=20, cex=.8, cex.axis=2, main='Stochastic')
plot(deterministic, pch=20, cex=.8, cex.axis=2, main='Deterministic')
#############################################
# Modelagem baseada em Processos Estocásticos
#############################################
#############################################
# Modelagem baseada em Sistemas Dinâmicas
#############################################
# Deterministic
par(mfrow=c(3,1))
# Estimar o time lag APENAS SOBRE A PARTE DETERMINISTIAS
ami = tseriesChaos::mutual(deterministic, lag.max=100)
v = diff(ami)
d = as.numeric(which(v > 0)[1])
cat("Time lag = ", d, "\n")
# Estimando Takens embedding dimension
fnn = tseriesChaos::false.nearest(deterministic, m=10, d=d, t=10)
plot(fnn)
print(fnn)
m = as.numeric(which.min(fnn[1,]))
# Aplicação do Takens embedding theorem
dataset = tseriesChaos::embedd(deterministic, m=m, d=d)
# Modelagem e Predição
labelId = ncol(dataset)
X = matrix(dataset[,1:(labelId-1)], ncol=labelId-1)
Y = dataset[,labelId]
buffer = dataset
for (new.row in (nrow(buffer) + 1) :(nrow(buffer) + 300)) {
x = as.numeric(buffer[new.row - d, 2:ncol(buffer)])
y = dwnn(query=x, X=X, Y=Y, sigma=sigma)
buffer = rbind(buffer, c(x,y))
}
plot(Y, t='l', cex.lab=2, cex=2, cex.axis=2,
xlab="Tempo", ylab='Valores', xlim=c(1, nrow(buffer)))
points(buffer[, labelId], col=2, pch=20, cex=.2)
return (valores)
}
x = test.seno_com_muito_ruido.decomposicao.espaco_fase()
No ponto em que há uma transição (comparação entre a 4 e 5 decomposição). Neste ponto, é onde começa a parte determinística
Como modelar o componente estocástico??
test.seno_com_muito_ruido.decomposicao.espaco_fase <- function(sigma=0.5, sd=0.5, mi.epsilon=0.05) {
tempo = seq(0,5, length=1000)
valores = sin(2*pi*tempo) + rnorm(mean=0, sd=sd, n=1000)
# Componente estocástico
require(EMD)
decomposicao = emd(valores)
phases.vector = list()
for (i in 1:decomposicao$nimf) {
coeff = fft(decomposicao$imf[,i])
phases = atan(Im(coeff) / Re(coeff))
phases.vector[[i]] = phases
}
coeff = fft(decomposicao$residue)
phases = atan(Im(coeff) / Re(coeff))
phases.vector[[i+1]] = phases
require(c3net)
mutual.information = rep(0, length(phases.vector)-1)
for (i in 1:(length(phases.vector) - 1)) {
mutual.information[i] = makemim(rbind(phases.vector[[i]], phases.vector[[i+1]]))[1,2]
}
plot(mutual.information, pch=20, cex=.8)
det.start = which(mutual.information > mi.epsilon)[1]
stochastic = rowSums(decomposicao$imf[,1:(det.start-1)])
# Componente determinístico
deterministic = rowSums(decomposicao$imf[,det.start:decomposicao$nimf]) + decomposicao$residue
par(mfrow=c(3,1))
plot(valores, pch=20, cex=.8, cex.axis=2, main='Serie temporal original')
plot(stochastic, pch=20, cex=.8, cex.axis=2, main='Stochastic')
plot(deterministic, pch=20, cex=.8, cex.axis=2, main='Deterministic')
#############################################
# Modelagem baseada em Processos Estocásticos
#############################################
# Stochastic
require(forecast)
model = forecast::auto.arima(stochastic)
stochastic.pred = as.numeric(predict(model, n.head=300)$pred)
#############################################
# Modelagem baseada em Sistemas Dinâmicas
#############################################
# Deterministic
par(mfrow=c(3,1))
# Estimar o time lag APENAS SOBRE A PARTE DETERMINISTIAS
ami = tseriesChaos::mutual(deterministic, lag.max=100)
v = diff(ami)
d = as.numeric(which(v > 0)[1])
cat("Time lag = ", d, "\n")
# Estimando Takens embedding dimension
fnn = tseriesChaos::false.nearest(deterministic, m=10, d=d, t=10)
plot(fnn)
print(fnn)
m = as.numeric(which.min(fnn[1,]))
# Aplicação do Takens embedding theorem
dataset = tseriesChaos::embedd(deterministic, m=m, d=d)
# Modelagem e Predição
labelId = ncol(dataset)
X = matrix(dataset[,1:(labelId-1)], ncol=labelId-1)
Y = dataset[,labelId]
buffer = dataset
for (new.row in (nrow(buffer) + 1) :(nrow(buffer) + 300)) {
x = as.numeric(buffer[new.row - d, 2:ncol(buffer)])
y = dwnn(query=x, X=X, Y=Y, sigma=sigma)
buffer = rbind(buffer, c(x,y))
}
deterministic.pred = buffer[(nrow(X)+1):nrow(buffer), labelId]
plot(Y, t='l', cex.lab=2, cex=2, cex.axis=2,
xlab="Tempo", ylab='Valores', xlim=c(1, nrow(buffer)))
stochastic.pred = c(rep(NA, nrow(X)), stochastic.pred)
points(stochastic.pred, col=2, pch=20, cex=.8)
plot(Y, t='l', cex.lab=2, cex=2, cex.axis=2,
xlab="Tempo", ylab='Valores', xlim=c(1, nrow(buffer)))
deterministic.pred = c(rep(NA, nrow(X)), deterministic.pred)
points(deterministic.pred, col=2, pch=20, cex=.8)
return (valores)
}
stochastic = test.seno_com_muito_ruido.decomposicao.espaco_fase()
Assim como FNN/AMI ajuda a encontrar Takens embedding theorem...
Para processos estocásticos: ACF/PACF, com o modelo ARIMA proposto por Box eJenkins
Ainda considera que existam um número de componentes determinísticas (com a letra $a$, AUTO REGRESSIVE AR - dependência com o passado), e um número de componentes estocásticas (com a letra $b$, MOVEL AVERAGE - processos estocásticos). Tomando a definição do manual da função 'arima' (I - INTEGRATES, um nível de dependência temporal entre as duas variáveis) inata de R:
$$ X[t] = a[1]X[t-1] + \ldots + a[p]X[t-p] + e[t] + b[1]e[t-1] + \ldots + b[q]e[t-q] $$install.packages('forecast')
require(forecast)
forecast::auto.arima(stochastic)
Com uma AR não nulo, existe uma relação determinística. Já com MA, existem três componentes com influencia, sendo a imediatamente anterior da medida (ma1), a mais relevante
Aplicação NOAA¶
paris = read.csv("Estudos/NOAAparisData/NOAAparis.csv")
row.ids = which(paris$STATION == "FRM00007149")
orly = paris[row.ids,]
colnames(orly)
plot(orly$TAVG)
sum(is.na(orly$TAVG[12000:25000]))
Significa que 46 valores neste intervalo não estão preenchidos
which(is.na(orly$TAVG[12000:25000]))
dwnn <- function(query, X, Y, sigma) {
E = apply(X, 1, function(row) { sqrt(sum((row - query)^2)) })
# Construção de gaussiana centrada na query
weight = exp(-E^2 / (2*sigma^2))
# Multiplicando ponto a ponto os pesos pela saída de cada elemento
return (weight %*% Y / sum(weight))
}
noaa <- function(filename='Estudos/NOAAparisData/NOAAparis.csv', station='FRM00007149', field='TAVG',
data.gap=0.01, k=25, sigma=2, mi.epsilon=0.001) {
dataset = read.csv(filename)
row.ids = which(dataset$STATION == station)
values = dataset[row.ids, field]
filled.ids = which(!is.na(values))
region = which(diff(filled.ids) > data.gap*length(values))
series = values[filled.ids[region[1]:region[2]]]
#plot(series)
# Preencimento devalores faltantes
# ---- séries temporais: interpolação da vizinhança
for (i in which(is.na(series))) {
neighbors = seq(i-k, i+k)
neighbors = neighbors[neighbors >= 1 & neighbors <= length(series)]
neighbors = neighbors[which(!is.naseries[neighbors])]
series[i] == mean(series[neighbors])
}
# points(series, col=2)
# Componente estocástico
require(EMD)
decomposicao = emd(series)
phases.vector = list()
for (i in 1:decomposicao$nimf) {
coeff = fft(decomposicao$imf[,i])
phases = atan(Im(coeff) / Re(coeff))
phases.vector[[i]] = phases
}
coeff = fft(decomposicao$residue)
phases = atan(Im(coeff) / Re(coeff))
phases.vector[[i+1]] = phases
require(c3net)
mutual.information = rep(0, length(phases.vector)-1)
for (i in 1:(length(phases.vector) - 1)) {
mutual.information[i] = makemim(rbind(phases.vector[[i]], phases.vector[[i+1]]))[1,2]
}
plot(mutual.information, pch=20, cex=.8)
det.start = which(mutual.information > mi.epsilon)[1]
stochastic = rowSums(decomposicao$imf[,1:(det.start-1)])
# Componente determinístico
deterministic = rowSums(decomposicao$imf[,det.start:decomposicao$nimf]) + decomposicao$residue
par(mfrow=c(3,1))
plot(series, pch=20, cex=.8, cex.axis=2, main='Serie temporal original')
plot(stochastic, pch=20, cex=.8, cex.axis=2, main='Stochastic')
plot(deterministic, pch=20, cex=.8, cex.axis=2, main='Deterministic')
#############################################
# Modelagem baseada em Processos Estocásticos
#############################################
# Stochastic
require(forecast)
model = forecast::auto.arima(stochastic)
stochastic.pred = as.numeric(predict(model, n.head=300)$pred)
#############################################
# Modelagem baseada em Sistemas Dinâmicas
#############################################
# Deterministic
par(mfrow=c(3,1))
# Estimar o time lag APENAS SOBRE A PARTE DETERMINISTIAS
ami = tseriesChaos::mutual(deterministic, lag.max=100)
v = diff(ami)
d = as.numeric(which(v > 0)[1])
cat("Time lag = ", d, "\n")
# Estimando Takens embedding dimension
fnn = tseriesChaos::false.nearest(deterministic, m=10, d=d, t=10)
plot(fnn)
print(fnn)
m = as.numeric(which.min(fnn[1,]))
# Aplicação do Takens embedding theorem
dataset = tseriesChaos::embedd(deterministic, m=m, d=d)
# Modelagem e Predição
labelId = ncol(dataset)
X = matrix(dataset[,1:(labelId-1)], ncol=labelId-1)
Y = dataset[,labelId]
buffer = dataset
for (new.row in (nrow(buffer) + 1) :(nrow(buffer) + 300)) {
x = as.numeric(buffer[new.row - d, 2:ncol(buffer)])
y = dwnn(query=x, X=X, Y=Y, sigma=sigma)
buffer = rbind(buffer, c(x,y))
}
deterministic.pred = buffer[(nrow(X)+1):nrow(buffer), labelId]
plot(stochastic, t='l', cex.lab=2, cex=2, cex.axis=2,
xlab="Tempo", ylab='Series', xlim=c(1, nrow(buffer)))
stochastic.pred = c(rep(NA, length(stochastic)), stochastic.pred)
points(stochastic.pred, col=2, pch=20, cex=.8)
plot(deterministic, t='l', cex.lab=2, cex=2, cex.axis=2,
xlab="Tempo", ylab='Series', xlim=c(1, nrow(buffer)))
deterministic.pred = c(rep(NA, length(deterministic)), deterministic.pred)
points(deterministic.pred, col=2, pch=20, cex=.8)
ret = list()
ret$series = series
ret$stochastic = stochastic
ret$deterministic = deterministic
ret$stochastic.pred = stochastic.pred
ret$deterministic.pred = deterministic.pred
ret$prediction = stochastic.pred + deterministic.pred
return (ret)
}
x = noaa()
Aula 07¶
Retomada dos algoritmos de IBL
knn <- function(query, k, X, Y) {
E = apply(X, 1, function(row) { sqrt(sum((row - query)^2)) })
row.ids = sort.list(E, dec=F)[1:k]
classes = unique(Y)
count = rep(0, length(classes))
i = 1
for (class in classes) {
count[i] = sum(class == Y[row.ids])
i = i+1
}
ret = list()
ret$classes = classes
ret$count = count
ret$max.prob.class = classes[which.max(count)]
return (ret)
}
dwnn <- function(query, X, Y, sigma) {
E = apply(X, 1, function(row) { sqrt(sum((row - query)^2)) })
weight = exp(-E^2 / (2*sigma^2))
return (weight %*% Y / sum(weight))
}
seno <- function(sigma=0.5, mi.epsilon=0.05) {
tempo = seq(0, 5, length=1000)
valores = sin(2*pi*tempo) + rnorm(mean=0, sd=0.15, n=1000)
# Decomposição: Applyind Empirical Mode Decomposition and mutual information
# to separate stochastic and deterministic influences embedded in signals
require(EMD)
decomposicao = emd(valores)
# decomposicao agora passa a carregar outras senóides que compõe a senóide principal
par(mfrow=c(3, ceiling((decomposicao$nimf + 2) / 3)))
plot(valores, pch=20, cex=.8, main='Série original', cex.axis=2, cex.main=2)
for (i in 1:decomposicao$nimf) {
plot(decomposicao$imf[,i], pch=20, cex=.8, main=paste("IMF", i), cex.axis=2, cex.main=2)
}
plot(decomposicao$residue, pch=20, cex=.8, main='Residuo', cex.axis=2, cex.main=2)
par(mfrow=c(3, ceiling((decomposicao$nimf + 2) / 3)))
plot(valores, pch=20, cex=.8, main='Série original', cex.axis=2, cex.main=2)
for (i in 1:decomposicao$nimf) {
coeff = fft(decomposicao$imf[,i])
plot(coeff, asp=1, pch=20, cex=.8, main=paste("IMF", i), cex.axis=2, cex.main=2)
}
coeff = fft(decomposicao$residue)
plot(coeff, asp=1, pch=20, cex=.8, main='Residuo', cex.axis=2, cex.main=2)
phases.vector = list()
par(mfrow=c(3, ceiling((decomposicao$nimf + 2) / 3)))
plot(valores, pch=20, cex=.8, main='Série original', cex.axis=2, cex.main=2)
for (i in 1:decomposicao$nimf) {
coeff = fft(decomposicao$imf[,i])
phases = atan(Im(coeff) / Re(coeff))
plot(phases, pch=20, cex=.8, main=paste("IMF", i), cex.axis=2, cex.main=2)
phases.vector[[i]] = phases
}
coeff = fft(decomposicao$residue)
phases = atan(Im(coeff) / Re(coeff))
plot(phases, pch=20, cex=.8, main='Residuo', cex.axis=2, cex.main=2)
phases.vector[[i+1]] = phases
# A congruência está relacionada com o casamento dos angulos dos vetores que compoe o espaço
# complexo de fourier, o que mostra que as senóides começam a ter comportamento mais suave de transição
require(c3net)
mutual.information = rep(0, length(phases.vector)-1)
for (i in 1:(length(phases.vector)-1)) {
mutual.information[i] = makemim(rbind(phases.vector[[i]], phases.vector[[i+1]]))[1,2]
}
# Makemim Compara as fases de cada decomposicao das imf. Se as fases contiverem grande quantidade
# de informação mútua, significa que começamos a ter congruência de fase, e começam a ter comportamento
# mais determinístico
# Ou também, qual o momento que as Intrinsic Mode Functions começam a ter similaridades relevantes entre si
par(mfrow=c(1,1))
plot(mutual.information, cex=.8, cex.axis=2, cex.main=1, main='Mutual Information', pch=20)
det.start = which(mutual.information > mi.epsilon)[1]
# Componente estocástico
stochastic = rowSums(decomposicao$imf[,1:(det.start-1)]) # Soma apenas das IMF mais estocásticas
#Componente determinístico
deterministic = rowSums(decomposicao$imf[,det.start:decomposicao$nimf]) + decomposicao$residue
par(mfrow=c(3,1))
plot(valores, pch=20, cex=.8, cex.axis=3, main='Serie temporal original')
plot(stochastic, pch=20, cex=.8, cex.axis=3, main='Stochastic')
plot(deterministic, pch=20, cex=.8, cex.axis=3, main='Deterministic')
ret = list()
ret$series = valores
ret$stochastic = stochastic
ret$deterministic = deterministic
return(ret)
}
result = seno()
Problema: Como modelar o cenário estocástico
Ferramentas: ACF (Auto-correlation Function), PACF(Partial Auto-Correlation Function)
(com a letra $a$, AUTO REGRESSIVE AR - dependência com o passado), e um número de componentes estocásticas (com a letra $b$, MOVEL AVERAGE - processos estocásticos). Tomando a definição do manual da função 'arima' (I - INTEGRATES, um nível de dependência temporal entre as duas variáveis) inata de R:
$$ X[t] = a[1]X[t-1] + \ldots + a[p]X[t-p] + e[t] + b[1]e[t-1] + \ldots + b[q]e[t-q] $$$p$: ordem do modelo auto-regressivo (AR)
$e$: erros associados a cada instante de tempo no passado
$q$: ordem do modelo de média móvel (MA)
O modelo ARIMA utiliza do modelo ARMA, mas realiza antes, uma decomposição da série temporal por meio de uma análise de estacionariedade
plot(result$stochastic, pch=20, cex=.8)
Este resultado (componente estocástico) é do tipo ESTACIONÁRIO, já que não existe uma tendência de crescimento ou decrescimento. Ou também os momentos estatísticos permanecem estáveis ao longo do tempo
mean(result$stochastic[1:500])
mean(result$stochastic[500:1000])
sqrt(sum((mean(result$stochastic[1:500]) - mean(result$stochastic[501:1000]))^2))
^ Calcula a diferença absoluta (distância euclidiana). Conforme aproxima de zero, temos que o primeiro momento (média) é estável ao longo das observações
Agora compararemos os quartos da série
sqrt(sum((mean(result$stochastic[1:250]) - mean(result$stochastic[251:500]))^2))
sqrt(sum((mean(result$stochastic[251:500]) - mean(result$stochastic[501:750]))^2))
sqrt(sum((mean(result$stochastic[500:750]) - mean(result$stochastic[751:1000]))^2))
Todos são próximos a zero!!!!
install.packages("moments")
estacionariedade <- function(series, threshold=250, order=1) {
require(moments)
# func: representa o momento a ser calculado (média, variância, curtose, skewness...)
len = length(series)
# A subtração das duas separações centraliza os momentos calculados. Por exemplo centraliza a média
# e podemos comparar deacordo com a proximidade de zero
div = sqrt((moment(series[1:floor(len/2)], order=order) - moment(series[(floor(len/2)+1):len], order=order))^2)
if (len > threshold) {
div = c(div, estacionariedade(series[1:floor(len/2)], order=order) +
estacionariedade(series[(floor(len/2)+1):len]), order=order)
}
return (div)
}
estacionariedade(result$stochastic)
^ a média varia próximo a zero
estacionariedade(result$stochastic, order=2)
^ a variância também
estacionariedade(result$stochastic, order=3)
novo.estocastico = result$stochastic + seq(0, 2, length=1000)
par(mfrow=c(2,1))
plot(result$stochastic, pch=20, cex=.8)
plot(novo.estocastico, pch=20, cex=.8)
O novo estocástico apresenta incremento ao longo do tempo, e portanto tem variação da média ao longo do tempo, e portanto não é mais estável
estacionariedade(novo.estocastico, order = 1)
Agora a média já varia bastante da proximidado do zero, o que indica que o componente não é estacionário
estacionariedade(novo.estocastico, order = 2)
estacionariedade(novo.estocastico, order = 3)
estacionariedade(novo.estocastico, order = 4)
Sistemas não estacionários possuem ordem de disparidade, ou nível. A estacionariedade teórica é nula, porém o nível
Os momentos estatísticos:: correlações entre eixos de um espaço¶
Conforme aumentamos o momento, podemos avaliar perturbações em um dado sistema ao longo do tempo por conta de que: caso o momento seja menor que 1, o aumento da ordem tende a reduzir seu valor. Caso contrário, o momento cresce
diff.primeira.ordem = diff(novo.estocastico)
estacionariedade(diff.primeira.ordem, order = 1)
diff.primeira.ordem = diff(novo.estocastico)
estacionariedade(diff.primeira.ordem, order = 2)
Quando calculamos a diferença ponto a ponto do atual com o anterior, removemos a tendência. Com isso podemos fazer a modelagem utilizando ARMA, com modelagem linear
series = arima.sim(n=1000, list(ar=c(0.8897))) # com ar, temos dependência apenas com passado próximo
#series[t] = a*series[t-1] + b
plot(series)
Com isso, esperamos que o espaço de fase entre $series[t]$ $vs$ $series[t-1]$ seja linear
require(tseriesChaos)
plot(embedd(series, m=2, d=1), pch=20, cex=.8)
Que de fato, obtemos
# A medida que aumentamos o lag com que plotamos, os pontos se dispersam devido ao acúmulo de erro
# Quando fizermos a projeção dos pontos sobre a regressão, para d=1 teremos menos erros
plot(embedd(series, m=2, d=1), pch=20, cex=.8)
points(embedd(series, m=2, d=2), pch=20, cex=.4, col=2)
points(embedd(series, m=2, d=3), pch=20, cex=.4, col=3)
points(embedd(series, m=2, d=8), pch=20, cex=.4, col=4)
require(stats)
dataset = embedd(series, m=2, d=1)
dataset = as.data.frame(embedd(series, m=2, d=1))
colnames(dataset) = c("Xt1", "Xt")
lm(formula=Xt~Xt1, data=dataset)
Xt1 sendo 0.875 representa o argumento com o qual construímos o modelo... o que faz sentido
series = arima.sim(n=1000, list(ar=c(0.8897, -0.25))) # Agora a dependência vai até a um lag com 2 passos
dataset = embedd(series, m=3, d=1)
dataset = as.data.frame(embedd(series, m=3, d=1))
colnames(dataset) = c("Xtm2", "Xtm1", "Xt")
lm(Xt ~., data=dataset)
install.packages('sourcetools')
series = arima.sim(n=1000, list(order=c(1,1,0), ar=0.7))
# a ordem do modelo será: c(1,1,0) - a ordem do AR é 1
# - a ordem da integração é 1
# - a ordem do MA é 0
Criamos assim uma série temporal na forma
Como a integração vale 1, a série é não estacionária e depende de uma dependência de primeira ordem
$$ series[t] = 0.7 * series[t-1] + b $$plot(series)
cat(estacionariedade(series, order=1))
cat(estacionariedade(series, order=2))
cat(estacionariedade(series, order=3))
cat(estacionariedade(series, order=4))
A ordem é MUITO GRANDE. O que faz da serie não estacionária
diff.primeira.ordem = diff(series)
par(mfrow=c(2,1))
plot(series)
plot(diff.primeira.ordem)
print(estacionariedade(diff.primeira.ordem, order=1))
print(estacionariedade(diff.primeira.ordem, order=2))
print(estacionariedade(diff.primeira.ordem, order=3))
print(estacionariedade(diff.primeira.ordem, order=4))
Agora podemos modelar a série utilizando ARMA
test <- function() {
series = arima.sim(list(order = c(1,2,0), ar=0.7), n = 1000)
# Agora temos integração de ordem 2, o que faz necessário a diferença de segunda ordem
cat("Calculando estacionariedades:")
print(estacionariedade(series, order=1))
# print(estacionariedade(series, order=2))
# print(estacionariedade(series, order=3))
# print(estacionariedade(series, order=4))
diff.primeira.ordem = diff(series)
cat("\nCalculando diferença de ordem 1:")
print(estacionariedade(diff.primeira.ordem, order=1))
# print(estacionariedade(diff.primeira.ordem, order=2))
# print(estacionariedade(diff.primeira.ordem, order=3))
# print(estacionariedade(diff.primeira.ordem, order=4))
diff.segunda.ordem = diff(diff.primeira.ordem)
cat("\nCalculando diferença de ordem 2:")
print(estacionariedade(diff.segunda.ordem, order=1))
# print(estacionariedade(diff.segunda.ordem, order=2))
# print(estacionariedade(diff.segunda.ordem, order=3))
# print(estacionariedade(diff.segunda.ordem, order=4))
par(mfrow=c(3,1))
plot(series)
plot(diff.primeira.ordem)
plot(diff.segunda.ordem)
return(diff.segunda.ordem)
}
diff2 = test()
Com estas etapas, descobrimos quais as dependências que precisamos remover antes da modelagem
Agora podemos calcular PACF:
pacf(diff2)
Isso significa que $lag = 1$ atraso no passado é o lag que mais influencia no evento atual dos 30 primeiros componentes analisados
test <- function() {
series = arima.sim(list(order = c(1,2,0), ar=0.7), n = 1000)
# Agora temos integração de ordem 2, o que faz necessário a diferença de segunda ordem
# cat("Calculando estacionariedades:")
# print(estacionariedade(series, order=1))
# print(estacionariedade(series, order=2))
# print(estacionariedade(series, order=3))
# print(estacionariedade(series, order=4))
diff.primeira.ordem = diff(series)
# cat("\nCalculando diferença de ordem 1:")
# print(estacionariedade(diff.primeira.ordem, order=1))
# print(estacionariedade(diff.primeira.ordem, order=2))
# print(estacionariedade(diff.primeira.ordem, order=3))
# print(estacionariedade(diff.primeira.ordem, order=4))
diff.segunda.ordem = diff(diff.primeira.ordem)
# cat("\nCalculando diferença de ordem 2:")
# print(estacionariedade(diff.segunda.ordem, order=1))
# print(estacionariedade(diff.segunda.ordem, order=2))
# print(estacionariedade(diff.segunda.ordem, order=3))
# print(estacionariedade(diff.segunda.ordem, order=4))
par(mfrow=c(3,1))
plot(series)
plot(diff.primeira.ordem)
plot(diff.segunda.ordem)
#pacf(diff2, lag.max=100)
return(diff.segunda.ordem)
}
diff2 = test()
pacf(diff2, lag.max=100)
acf(diff2, lag.max=100)
A ACF mostra que o erro com $lag = 5$ anteriores, ou seja, a ordem de MA é $5$
No modelo simulado, impomos um lag zero. A ACF é um estimador e existe um erro associado à própria ACF
Precisamos então testar uma lista de modelos candidatos: AR = 1; MA = 5 Logo:
(1,2,0); (1,2,1); (1,2,2); (1,2,3); (1,2,4); (1,2,5)
(0,2,0); (0,2,1); (0,2,2); (0,2,3); (0,2,4); (0,2,5)
arima(series, order=c(1,2,0))
Com estes parâmetros, o coeficiente encontrado pela função é de $-0.21$. Bem distante do esperado para $0.7$ (usado para construir o modelo)
arima(series, order=c(1,0,0))
arima(series, order=c(1,1,0))
arima(series, order=c(1,2,1))
arima(series, order=c(1,2,0))
aic: Akaike information index
sigma^2: erro associado ao modelo que queremos
Iremos escolher (1,2,1)
series = arima.sim(list(order = c(1,2,0), ar=0.7), n = 1000)
model = arima(series, order=c(1,2,1))
plot(c(series, rep(NA, 100)), pch=20, cex=.8)
points(c(rep(NA, length(series)), predict(model, n.ahead=100)$pred), col=2, pch=20, cex=.8)
test <- function() {
series = arima.sim(list(order = c(1,2,0), ar=0.7), n = 1000)
# Agora temos integração de ordem 2, o que faz necessário a diferença de segunda ordem
# cat("Calculando estacionariedades:")
# print(estacionariedade(series, order=1))
# print(estacionariedade(series, order=2))
# print(estacionariedade(series, order=3))
# print(estacionariedade(series, order=4))
diff.primeira.ordem = diff(series)
# cat("\nCalculando diferença de ordem 1:")
# print(estacionariedade(diff.primeira.ordem, order=1))
# print(estacionariedade(diff.primeira.ordem, order=2))
# print(estacionariedade(diff.primeira.ordem, order=3))
# print(estacionariedade(diff.primeira.ordem, order=4))
diff.segunda.ordem = diff(diff.primeira.ordem)
# cat("\nCalculando diferença de ordem 2:")
# print(estacionariedade(diff.segunda.ordem, order=1))
# print(estacionariedade(diff.segunda.ordem, order=2))
# print(estacionariedade(diff.segunda.ordem, order=3))
# print(estacionariedade(diff.segunda.ordem, order=4))
plot(series)
plot(diff.primeira.ordem)
plot(diff.segunda.ordem)
model = list()
model[[1]] = arima(x=series, order=c(0,2,0))
model[[2]] = arima(x=series, order=c(0,2,1))
model[[3]] = arima(x=series, order=c(0,2,2))
model[[4]] = arima(x=series, order=c(0,2,3))
model[[5]] = arima(x=series, order=c(0,2,4))
model[[6]] = arima(x=series, order=c(0,2,5))
model[[7]] = arima(x=series, order=c(1,2,0))
model[[8]] = arima(x=series, order=c(1,2,1))
model[[9]] = arima(x=series, order=c(1,2,2))
model[[10]] = arima(x=series, order=c(1,2,3))
model[[11]] = arima(x=series, order=c(1,2,4))
model[[12]] = arima(x=series, order=c(1,2,5))
model.id = -1
akaike.information = Inf
for (i in 1:length(model)) {
if (model[[i]]$aic < akaike.information) {
model.id = i
akaike.information = model[[i]]$aic
}
}
plot(pacf(diff.segunda.ordem, lag.max=100))
plot(acf(diff.segunda.ordem, lag.max=100))
cat("Selected model: ", model.id)
selected.model = model[[model.id]]
plot(c(series, rep(NA, 100)), pch=20, cex=.8)
points(c(rep(NA, length(series)), predict(selected.model, n.ahead=100)$pred), col=2, pch=20, cex=.8)
return(diff.segunda.ordem)
}
result = test()
Aula 08¶
OVERVIEW
- Dado dataset:
- Pegar um dado atributo
- Analisar dependência temporal do atributo (ACF)
- Caso dependente: Decomposição (EMD)
- Componente Determinístico: Modelo de sistemas dinâmicos
- AMI (Time lag)
- FNN (embedding)
- Takens' embedding theorem (traz independência aos atributos)
- Algoritmo de regressão (Reg. logística, DWNN, Support Vector regression, multilayer perceptron, RBF)
- Componente Determinístico: Modelo de sistemas dinâmicos
2. Componente Estocástico: Modelo de processos estocásticos
1. Análise de estacionariedade
2. PACF (para modelo AR - Auto-regressive)
3. ACF (para modelo de MA - Movel Average) (Reconstrução de espaço, que traz independência)
4. ARIMA (Fitting ou modelagem)
3. Caso independente (só é independente se ):
1. Aprendizado de máquina: Statistical Learning Theory (só tem validade neste caso onde as variáveis são independentes)dataset = read.table("Estudos/WineData/wine.data", sep=',')
acf(dataset[,2])
De acordo com esta análise, a ACF mostra que temos dependência entre os valores de um único atributo do dataset wine, o que não é verdade
Como não existem muitos dados, a ACF não traz muita precisão. Logo, saber se é dependente usando ACF só deve ser feito caso tenhamos grande quantidade de observações no conjunto de dados
Caso não existam muitos dados, podemos fazer
- segmentação: subconjunto treinamento/teste K vezes
- Avaliar a distribuiçãode acurácias
- Caso tenham uma distribuição de acuráci em torno da média, temos valores independentes
acf(sample(dataset[,2]))
Qualquer aleatorização quebra as dependências temporais.
dataset = read.table("Estudos/CovertypeData/covertype_csv.csv", sep=',')
print(dataset[1,])
Não existe correlação entre cada um dos valores
acf(dataset[,'V1'], lag.max=500)
Mesmo com muitos dados, existe uma dependência entre elas pela ACF, apesar de não existir. ISso ocorre porque a elevação foi ordenada
acf(sample(dataset[, 'V1']), lag.max=500)
Como antes, aleatorização quebra a dependência temporal
A Statistical Learning Theory é nossa garantia de que o aprendizado de máquina pode acontecer. Temos então o Modelo-f
Construimos o modelo f usando uma subamostra do dataset. O Poder de generalização do nosso modelo é dado pela diferença de resultado do modelo-f com o conjunto de treino com 75% das amostras em relação ao resultado do modelo-f com o conjunto de teste com 25% das amostras
Quão mais próximo de zero, maior a generalização. Para que SLT garanta aprendizado:
- Deve ser possível generalizar (R_emp1(f) - R_emp2(f) -> 0) -- R_emp: Erro da amostra empírica 1
- O risco na amostra (1 - R__emp1) "complementar de R_emp1" de treinamento deve tender a zero conforme o tamanho
Logo, tanto o poder de generalização quanto os resultados do modelo-f devem tender a zero conforme a amostra aumente
E somente funciona SE E SOMENTE SE:
- Os exemplos forem independentes entre si
- ...
install.packages('randomForest')
dwnn <- function(query, X, Y, sigma) {
E = apply(X, 1, function(row) {sqrt(sum((row - query)^2)) })
weight = exp(-E^2 / (2*sigma^2))
return (weight %*% Y / sum(weight))
}
x = sin(2*pi*seq(0,9,len=1000))
x_new = sample(x)
par(mfrow=c(2,2))
plot(x)
acf(x)
plot(x_new)
acf(x_new)
x = sin(2*pi*seq(0,9,len=1000))
x_new = sample(x)
x.phase.space = cbind(x[1:999], x[2:1000])
x_new.phase.space = cbind(x_new[1:999], x_new[2:1000])
x.buffer = x.phase.space[1:500,]
x_new.buffer = x_new.phase.space[1:500,]
sigma = 0.0015
# Modelagem e predição
for (i in (nrow(x.buffer) + 1):nrow(x.phase.space)) {
x = x.buffer[i-1, 2:ncol(x.buffer)]
y = dwnn(query=x, X=matrix(x.buffer[,1], ncol=1), Y=matrix(x.buffer[,2], ncol=1), sigma=sigma)
x.buffer = rbind(x.buffer, c(x,y))
}
for (i in (nrow(x_new.buffer) + 1):nrow(x_new.phase.space)) {
x = x_new.buffer[i-1, 2:ncol(x_new.buffer)]
y = dwnn(query=x, X=matrix(x_new.buffer[,1], ncol=1), Y=matrix(x_new.buffer[,2], ncol=1), sigma=sigma)
x_new.buffer = rbind(x_new.buffer, c(x,y))
}
par(mfrow=c(2,1))
plot(x.phase.space[,2], pch=20)
lines(x.buffer[,2], col=2, pch=20)
plot(x_new.phase.space[,2], pch=20)
lines(x_new.buffer[,2], col=2, pch=20)
x = sin(2*pi*seq(0,9,len=1000))
x_new = sample(x)
x.phase.space = cbind(x[1:999], x[2:1000])
x_new.phase.space = cbind(x_new[1:999], x_new[2:1000])
x.buffer = x.phase.space[1:500,]
x_new.buffer = x_new.phase.space[1:500,]
x.Error = NULL
x_new.Error = NULL
for (sigma in seq(0.000015, 0.001, length=10)){
x.buffer = x.phase.space[1:500,]
x_new.buffer = x_new.phase.space[1:500,]
# Modelagem e predição
for (i in (nrow(x.buffer) + 1):nrow(x.phase.space)) {
x = x.buffer[i-1, 2:ncol(x.buffer)]
y = dwnn(query=x, X=matrix(x.buffer[,1], ncol=1), Y=matrix(x.buffer[,2], ncol=1), sigma=sigma)
x.buffer = rbind(x.buffer, c(x,y))
}
Error = sqrt(sum((x.phase.space[501:nrow(x.phase.space),2] - x.buffer[501:nrow(x.phase.space),2])^2))
x.Error = rbind(x.Error, c(sigma,Error))
for (i in (nrow(x_new.buffer) + 1):nrow(x_new.phase.space)) {
x = x_new.buffer[i-1, 2:ncol(x_new.buffer)]
y = dwnn(query=x, X=matrix(x_new.buffer[,1], ncol=1), Y=matrix(x_new.buffer[,2], ncol=1), sigma=sigma)
x_new.buffer = rbind(x_new.buffer, c(x,y))
}
Error = sqrt(sum((x_new.phase.space[501:nrow(x_new.phase.space),2] - x_new.buffer[501:nrow(x_new.phase.space),2])^2))
x_new.Error = rbind(x_new.Error, c(sigma,Error))
}
par(mfrow=c(1,1))
x.ids = which(!is.nan(x.Error[,2]))
x_new.ids = which(!is.nan(x_new.Error[,2]))
plot(x.Error, ylim=range(c(x.Error[x.ids,2], x_new.Error[x_new.ids,2])), t='l')
lines(x_new.Error, col=2)
Para a faiza de valores de sigma bastante pequenos, o erro associado ao modelo com o espaço completo (senóide) é bastante baixo (próximo a zero, em preto). Já para o modelo com amostras aleatórias, o erro é sempre alto (vermelho)
Isso significa que a ACF apresentar dependencia temporal, usar DWNN somente não teremos aprendizado bom. Precisamos usar de Decomposição e análise determinística e estocástica
Aula 10¶
install.packages("class")
Para um aprendizado suprevisionado, para qualquer modelo f induzzido a partir de um conjunto de exemplos a fim de f: $X -> Y$
tal que x1, ..., xn $\in X$ é o espaço de entrada e y1, ..., yn $\in Y$ é o espaço de saída
Como garantir o sucesso do aprendizado supervisionado??
Vapnik -> Teoria do aprendizado estatístico: Conceito de Generalização ajuda nas provas de aprendizado
[Remp(f) - R(f)] $\to 0$ se n $\to \infty$
Onde
Remp(f) = (1/n) sum(i, 1, n, Loss(xi, yi, f(xi)))(Risco empírico, calculado sobre uma amostra apenas)
0-1 Loss: 0-1 Loss(xi, yi, f(xi)) = 0, se f(xi) = yi; = 1, caso constrário
R(f) = E(Loss(X, Y, f(X)) (Valor esperado de Loss), aplicado sobre todas as possíveis entradas Instanciação do conceito de generalização
Seja um professor responsável pelo ensono de tal assunto perante uma Classe
É responsabilidade do professor criar um conjunto de exemplos representativos para o problema em questão:
- Esses exemplos correspondem à realidade ou são um subconjunto que não contempla todas as possíveis classes (classificação) ou espaços de saída (regressão)?
- Motivar os alunos a procurar por materiais de viéses distintos (codificar, ler livros ou tutoriais com características distintas e diversas)
É responsabilidade de cada aluno:
- Não decorar (overfitting) um material único. O overfitting é a representação perfeita de uma amostra única, então Remp(f) = 0. Nisso, R(f) = 0 apenas se o problema analisado é simplório. Caso um problema mais real/complexo, R(f) = max
Uma prova (instância):
- $|Remp(f) - R(f)| \to 0$, se $n \to \infty$
- Nota na prova (acurácia): 1 - ((1/n) sum(i, 1, n, 0-1 Loss(xi, yi, f(xi)))). Nota 0.9, Remp 0.1; Nota 0.7, Remp 0.3
- Quando colocado em problemas do mundoreal, estes alunos com essas notas serão (Generalização de Vapnik): Mundo Real (Remp': indica que o mundo real não é um modelo generalista completo. Existem viéses). Remp' = 0.15; R'emp = 0.3
- A generalização será |0.1 - 0.15| = 0.05; |0.3 - 0.3| = 0
O que é APRENDER para Vapnik?? (Convergências)
- Generalização $|Remp(f) - R(f)| \to 0$, $n \to \infty$ (Precisamos ter uma prova com infinitos elementos para o aluno tenha um conhecimento generalizado)
- $Remp(f) \to 0$, $n \to \infty$ (Precisamos ter um aprendizado, dentro da própria amostra, deve ser perfeita- overfitting
A generalização é apenas um critério de seleção de modelos (estimador)
se $n \to \infty$, com $Remp(f) \in [0,1]$, $R(f) \in [0,1]$
Sendo f o viés (bias) do algoritmo de aprendizado utilizado
$$ P(sup_{f \in F} |Remp(f) - R(f)| > \epsilon) \leq 2 P(sup_{f \in F} |Remp(f) - Remp^*(f)| > \epsilon) \leq 2 m e^{-n\epsilon^2 / 4}) $$Onde m: número de funções contida no viés do algoritmo (numero de funções que o algoritmo é capaz de representar)
Hold-Out ---> Amostra 1 para treino | Amostra 2 para teste
- Se $|Remp(Amostra1) - Remp(Amostra2)| \to 0$; então generalizamos
- Se $Remp(Amostra1) \to 0$ ou $Remp(Amostra2) \to 0$; então aprendemos
A1 (treino) | A2 (validação) | A3 (teste)
overfitting | subseleção | teste
Subsampling ---> A1 (n) reamostrados com reposição 68\% (treino) | A2 disjunto de A1 (teste)
Em A2 temos exemplos nunca vistos em treino
K-fold Cross Validation ---> K folds ($k \leq n$). Se $k = n$, temos Leave one out
fold 1 (1/n) fold 2 (1/n) fold 3 (1/n) ... fold K (1/n)
Resto de folds (treino) | Fold $i$ (teste)
Algoritmo $A_i$:
- Remp(fold1)(f), ... Remp(foldK)(f)
- mean(Remp(fold1)(f), ... Remp(foldK)(f))
- sd(Remp(fold1)(f), ... Remp(foldK)(f))
Escolher o algoritmo $A_i$ com menor mean e sd
Um outro algoritmo B:
sendo
$$ 2 P(sup_{f \in F} |Remp(f) - Remp^*(f)| > \epsilon) \leq 2 m e^{-n\epsilon^2 / 4}) $$Definimos $\delta = 2 m e^{-n\epsilon^2 / 4}$ tal que
$$ \epsilon = \sqrt{\frac{4}{n} log(2m) - log(\delta)} $$Quando erramos: $$ sup_{f \in F} |Remp(f) - Remp'(f)| > \sqrt{\frac{4}{n} log(2m) - log(\delta)} $$
Quando acertamos: $$ sup_{f \in F} |Remp(f) - Remp'(f)| \leq \sqrt{\frac{4}{n} log(2m) - log(\delta)} $$
Se calcularmos uma matriz de divergências absolutas entre todos os riscos empríricos, teremos uma distribuição de erros. Quando \epsilon se aproxima de zero, temos o K-fold cross validation
Como o conjunto iris tem 50 setosa, 50 versicolor, 50 virginica
com K=10, temos ~Leave one out (0.9 K fold)
install.packages("dismo")
require(dismo)
folds.ids = kfold(as.numeric(iris[,5]), k=10)
print(folds.ids)
iris[folds.ids==2, ]
representam o id de cada fold onde colocamos uma das 10 amostras. Cada fold não terá o mesmo número de setosa, virginica e versicolor (não segue a mesma distribuição dos dados contidos no conjunto)
Precisamos então criar folds estratificdos, i.e., contendo a mesma representatividade dasclasses existentes no nosso conjunto de dados original
folds = list()
for (class in unique(iris[,5])) {
ids = which(iris[,5] == class)
folds.ids = kfold(ids, k=10)
for (fold in 1:10) {
sel = which(fold == folds.ids)
if (length(folds) >= fold){
folds[[fold]] = rbind(folds[[fold]], iris[ids[sel],])
} else {
folds[[fold]] = iris[ids[sel],]
}
}
}
folds[[1]]
Agora temos o mesmo número de elementos
# Precisamos randomizar cada fold
for (fold in 1:10) {
folds[[fold]] = folds[[fold]][sample(1:nrow(folds[[fold]])),]
}
folds[[1]]
require(class)
Remp.folds = c()
for (fold in 1:10) {
X.train = NULL
Y.train = c()
for (f in setdiff(1:10, fold)) {
X.train = rbind(X.train, folds[[f]][,1:4])
Y.train = c(Y.train, folds[[f]][,5])
}
y = knn(train=X.train, test=folds[[fold]][,1:4], cl=Y.train, k=3)
acc = sum((as.numeric(y) == as.numeric(folds[[fold]][,5])) / nrow(folds[[fold]]))
Remp = 1 - acc
Remp.folds = c(Remp.folds, Remp)
}
Remp.folds
cat('Remp médio = ', mean(Remp.folds), '\n')
cat('Remp stdev = ', sd(Remp.folds), '\n')
Cada modelo apresenta um resultado diferente do outro. Eventualmente os exemplos do conjuto de dados apresentam algum grau de dependência entre si
Pode ser que o desbalanceamentode classes por fold afeta o aprendizado (Se tivermos 90% de virginica, os kfolds terá esse viés, e o aprendizado será baseado apenas nela. Afeta muito a medida de risco Remp)
install.packages('tm')
Aula 11¶
require(tm)
library(tm)
news = "Estudos/20newsData"
dataset = VCorpus(DirSource(news, recursive=TRUE))
dataset
DTM = DocumentTermMatrix(dataset)
Geramos uma matriz com cada termo presente nos textos calculando a frequência de cada um dos termos
naive <- function(query, X, Y, eps=1e-323) {
H = unique(Y)
prob = rep(1, length(H))
P_h = rep(0, length(H))
for (i in 1:length(H)){
P_h[i] = sum(H[i] == Y) / nrow(X)
ids = which(uery > 0)
mat = as.matrix(X[,ids])
row.ids = which(rownames(mat) == X$dimnames$Docs)
P_aj_h = sum(H[i] == Y[row.ids]) / sum(H[i] == Y)
prob[i] = prob[i] * P_h[i] + eps
}
ret = list()
ret$H = H
ret$prob = prob / sum(prob)
return (ret)
}
dataset = VCorpus(DirSource(news, mode='text', recursive=TRUE))
strsplit_space_tokenizer <- function(x) {
unlist(strsplit(as.character(x), "[[:space:]]+"))
}
ctrl = list(tokenize=strsplit_space_tokenizer,
removePunctuation = TRUE,
stopwords=stopwords('en'),
stemming = FALSE,
tolower = TRUE,
removeNumbers=TRUE,
wordLengths=c(2,20))
X = DocumentTermMatrix(dataset, control = ctrl)
file.paths = system(paste('find ', news, ' -type f', sep=''), intern=TRUE)
Y = rep('none', nrow(X))
for (i in 1:length(X$dimnames$Docs)) {
id = grep(X$dimnames$Docs[i], file.paths)
elements = strsplit(file.paths[id], '/')[[1]]
Y[i] = elements[length(elements)-1]
}
Aula 12¶